module mrk4
    implicit none
    real,parameter :: h=0.05,h2=h/2.0,h6=h/6.0,xlim(2)=(/0.0,0.8/),y0=1.0
    integer,parameter :: N=int((xlim(2)-xlim(1))/h)
    real,save :: x(N),y(N)
    integer,save :: sta
contains
    real function f(x,y)
        implicit none
        real :: x,y
        f = y**2.0*cos(x)
    end function
    
    subroutine init
        implicit none
        integer :: i
        forall(i=1:N)x(i)=xlim(1)+h*(i-1)
        y=0.0
        y(1)=y0
        sta=1
    end subroutine
    
    subroutine rk4(func)
        implicit none
        real :: k(4)
        real,external :: func
        k(1)=func(x(sta),y(sta))
        k(2)=func(x(sta)+h2,y(sta)+h2*k(1))
        k(3)=func(x(sta)+h2,y(sta)+h2*k(2))
        k(4)=func(x(sta)+h,y(sta)+h*k(3))
        y(sta+1)=y(sta)+h6*(sum(k)+sum(k(2:3)))
        sta=sta+1
    end subroutine
end module

program main
    use mrk4
    implicit none
    real :: y2(N)
    open(233,file='data.dat',form='formatted',access='sequential')
    call init
    y2=1.0/(1.0-sin(x))
    write(233,100)x(sta),y(sta),y2(sta)
    do while(sta<N)
        call rk4(f)
        write(233,100)x(sta),y(sta),y2(sta)
    end do
    close(233)
    write(*,*)sin(3.14159/2.0)
100 format(f5.2,',',f12.9,',',f12.9)
end program
